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ABSTRACT 

We develop a new perturbative framework for studying the r-modes of rotating super- 
fluid neutron stars. Our analysis accounts for the centrifugal deformation of the star, 
and considers the two-fluid dynamics at linear order in the perturbed velocities. Our 
main focus is on a simple model system where the total density profile is that of an 
n = 1 polytrope. We derive a partially analytic solution for the superfluid analogue of 
the classical r-mode. This solution is used to analyse the relevance of the vortex medi- 
ated mutual friction damping, confirming that this dissipation mechanism is unlikely 
to suppress the gravitational-wave driven instability in rapidly spinning superfiuid 
neutron stars. Our calculation of the superfluid r-modes is significantly simpler than 
previous approaches, because it decouples the r-mode from all other inertial modes of 
the system. This leads to the results being clearer, but it also means that we cannot 
comment on the relevance of potential avoided crossings (and associated "resonances" ) 
that may occur for particular parameter values. Our analysis of the mutual friction 
damping differs from previous studies in two important ways. Firstly, we incorporate 
realistic pairing gaps which means that the regions of superfluidity in the star's core 
vary with temperature. Secondly, we allow the mutual friction parameters to take the 
whole range of permissible values rather than focussing on a particular mechanism. 
Thus, we consider not only the weak drag regime, but also the strong drag regime 
where the fluid dynamics is significantly different. 



1 INTRODUCTION 



Oscillations of rapidly rotating neutron stars have attracted interest for a considerable time. During the last decade significant 
effort was aimed at understanding whether gravitational-wave emission sets the upper speed limit for pulsars, e.g. via the r- 
mode instability ( Andersson] 1998 1. This possibility is of particular interest since such unstable systems may radiate detectable 
gravitational waves. In particular, it was suggested that r-modes in rapidly rotating neutron stars in low-mass X-ray binaries 



(LMXBs) may lead to a persistent source of gravitational radiation (Bildsten 1998 Andersson et al. 1999b I that may be 
detected by advanced LIGO [see [Watts et al.| ( |2008[ ) for the most recent analysis of this problem]. It is of great importance 
to understand whether internal fluid dissipation allows the instability to develop in such systems or whether it suppresses the 
r-modes completely. 

There are, in fact, many mechanisms at work in a real neutron star which compete with the gravitational- wave driving of 
the r-mode. One can express the relative strength of these mechanisms in terms of timescales which will, in general, depend 
on a large number of parameters. For a basic instability analysis, the most important parameters are the rotation rate and the 
temperature of the star. The instability can only develop when the gravitational radiation growth time-scale is shorter than the 
damping timescales due to the various viscosity mechanisms. This deflnes a region in the spin-temperature parameter space 
where the r-mode instability is active. Several studies have been devoted to calculating the shear- and bulk viscosity timescales 
for more or less realistic neutron star models. These studies tend to agree that the notion of persistent gravitational-wave 



emission accreting neutron stars in LMXBs is quite robust [see Nayyar & Owen ( 2006 1 and references therein] 



The neutron stars in LMXBs are believed to be old cold recycled pulsars that have been spun up by accretion. These 
stars are expected to have superfluid cores and are thus expected to have rather different fluid dynamics and dissipation 
mechanisms. The simplest description of these systems is provided by a two-fluid model, where the superfluid neutrons are 
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distinguished from the proton-electron plasma (usually treated as a charge neutral fluid). The oscillations of neutron stars with 
such superfiuid cores have been studied by Lindblom & Mendell (19941 and Andersson & Comer (20011. The latter showed 
that there are no g-modes in such stars, but on the other hand one has a new class of "superfiuid" modes. In particular, there 
are two families of r-modes ( Prix et ar]|2004 1 . One of these resembles the ordinary r-modes of a barotropic star in the sense 
that the neutrons and protons "move together" . For the other set of modes the neutrons and protons are "counter-moving" . 
The relative motion of neutrons and protons in a superfiuid is, however, rapidly damped by mutual friction, a mechanism that 
is directly associated with the rotational neutron vortices. The most commonly considered mutual friction mechanism is the 
scattering of electrons off the magnetic fields associated with the superfiuid neutron vortices ( Alpar et al.|1984 Andersson et 
|al.''2006). Mutual friction is known to have a decisive impact on neutron star dynamics. In fact, Lindblom & Mendell (19951 



have argued that mutual friction completely suppresses the gravitational radiation driven f-mode instability [see Andersson et 
|al.| ( |2008| for a recent discussion]. It is thus important to understand if the r-mode instability suffers a similar fate. Previous 
work on the subject by [Lindblom fc Mendell ( 2000 \ and Lee & Yoshida ( 2003 1 presents a somewhat mixed picture. While 
the ordinary r-mode is very weakly damped for most values of the entrainment parameter, the damping becomes very strong 
for a small range of parameter values. Hence, there may be situations where the mutual friction is strong enough to stop the 
instability from growing. 

The purpose of this paper is to clarify the role of the mutual friction damping and understand whether it can stabilize the 
r-modes against the gravitational- wave instability. In order to do this we extend the formulation of Andersson et al. ( 2008 1 to 
second order in rotation to calculate the mutual damping timescale for r-modes. Although we do not expect our analysis to 
lead to results that change previous conclusions, there are good reasons to return to this problem. First of all, previous studies 
have focussed entirely on the weak drag regime for the mutual friction. Meanwhile, recent studies of superfiuid turbulence 
Andersson et al. 2007: Peralta et al.|2005 20061, neutron star free precession ( Glampedakis et al.|2008a|b I and pulsar glitches 



Glampedakis fc Andersson||2008 1 demonstrate that systems in the strong drag regime may have very dynamics. Since one 



can make convincing arguments for the strong drag regime being relevant ( [Ruderman et al.|1998||Link|2003[p006l[SedrakIan| 
fc Sedrakian||1995 1, we clearly need to understand the r-mode problem for this parameter range as well. Secondly, we want 



to improve the treatment of the critical temperature/density at which superfluidity comes into play. The density dependence 
of the various superfiuid pairing gaps translates into distinct regions of normal- and super-fluids that vary with the core 
temperature. So far, the best analysis in this respect is that of Lindblom & Mendell (20001 who consider a two-fluid core 
surrounded by a single fluid envelope. We will build on this by considering superfluid layers, the size of which are determined 
by a (qualitatively) realistic pairing gap. The size of the superfluid region affects not only the mutual friction but also the 
shear- and bulk viscosities, and it is relevant to demonstrate how this alters the r-mode instability window. Finally, we want 
to lay the foundation for more detailed work on exotic neutron star cores. It is well established that exotic cores, dominated 



by either hyperons or deconfined quarks, may be associated with a very strong bulk viscosity [see, for example, Nayyar & 



Owen (20061]. Although it has been suggested that this would lead to the suppression of the r-modes, it is clear that such a 



conclusion is premature. In reality one would expect superfluidity to play a role. If the hyperons are superfluid the reactions 



that produce the bulk viscosity are suppressed and the effect on the r-mode instability may not be that severe (Nayyar & 



Owen|2006 1. A similar argument applies to colour superconducting quarks ( Alford et al.|1999 l. However, it may be important 



to keep in mind that a superfluid system has extra dynamical degrees of freedom. In order to truly understand the dynamics 
of these exotic phases of matter we need to explore the multi-fluid nature of these systems. The present study paves the way 
for future work in this direction. 



2 THE TWO-FLUID EQUATIONS 
2.1 The unperturbed problem 

Our discussion is based on the usual two-fluid model for neutron star cores ([Prix 2004 Andersson & Comer 20061. That 



is, we consider two dynamical degrees of freedom, loosely speaking representing the superfluid neutrons (labeled n) and a 
charge-neutral conglomerate of the protons and electrons (labeled p). Assuming that the individual species are conserved, we 
have the standard conservation laws 

atp, + V.(pX) = (1) 

where the constituent index x may be either p or n. The equations of momentum balance can be written 

(dt + vi\/j)iv: + e^wD + V.(/ix + <&) + ExKxViW," = /r/px (2) 

where = — vf (y 7^ x), and p.^ = ^^/m^ represents the chemical potential (in the following we assume that rrip = m^). 
Moreover, $ represents the gravitational potential, and the parameter £x encodes the entrainment effect. The force on the 
right-hand side of ([2| can be used to represent other interactions, including dissipative terms ( Andersson fc Comer|20(36 1. We 



will focus on the vortex-mediated mutual friction force for a system that, in equilibrium, rotates uniformly. This means that 
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we consider a force of form ( Andersson at al.|2006 l 

Here, fl-' is the angular frequency of the neutron fluid (a hat represents a unit vector). 

One can express the mutual friction force in terms of a dimensionless "drag" parameter TZ such that ( [Andersson et al" 



(3) 



20061 



B = 



TZ 



and 



B' = 



7^2 



(4) 



In the standard picture, the mutual friction is due to the scattering of electrons off of an array of neutron vortices ( Alpar et al 



19841 which leads to 7?. ^ 1, thus placing the problem in the weak drag regime. This means that B ^ 1 and B' <^ B. Hence, 



all B' terms can be neglected. However, it is commonly thought ( Ruderman et al.|l998 Link|20d3 20061 that the strong drag 
limit, 7?. 3> 1, may apply. Intermediate values for the drag, TZ ~ 1 are also interesting. In particular since one would then 
have both B ~ 1 and B' ~ 1 (note that in the case TZ ^ 1 one still has S' « 1 but B <^ 1) . This leads to the presence of 
terms that may have significant effect on the dynamics of the system ( Glampedakis et al.|2008a|b Glampedakis & Andersson 
20081. At this point, our understanding of neutron star core physics is sufHciently rudimentary that we should avoid ruling 



out the various possibilities. Hence, we will consider the entire permissible range of values for the drag parameter. 

Let us consider a frame rotating with the star at fixed angular velocity Q^. The equations of motion then take the form: 

(5) 



where we have included the centrifugal term in the potential 

<l>j? = $ — 

The continuity equations maintain the form |T]) and the Poisson equation 



1 o2 2 . 2 a 

-il r sm o 



(6) 



(7) 



2.2 Perturbations 

To keep the problem tractable we will assume that the background configuration is such that the two fiuids rotate together 



[see Prix et al. (20041; Glampedakis & Andersson ( ^2008) for discussions of oscillations in systems that are not in co-rotation]. 
Perturbing the equations of motion and working in a frame rotating with we then have 



and 



at5px + V,(pxK) = 



(8) 



(9) 



To completely specify the perturbation problem, we need boundary conditions. At the centre of the star we simply require 
that all variables are regular. The surface of the star is somewhat more complex. In reality one does not expect the superfluid 
region to extend all the way to the surface. We shall thus assume that the superfluid is only present in a distinct region 
determined by the core temperature and the superfluid pairing gap. We will discuss the implementation of this idea once we 
have set up the relevant system of equations. 

From previous work on superfluid neutron star oscillations, e.g. Lindblom & Mendell ( 19941; Andersson & Comer (2001 1; 



Prix & Rieutord (20021; Andersson et al. (20081, we know that the problem has two "natural" degrees of freedom. One of 

(10) 



them represents the total mass flux. Introducing 

pSv^ = PnSvi + ppSv^p 

and combining the two Euler equations accordingly, we get 



Sp 







(11) 



where p — pn 



p 

Pp and the pressure is obtained from 

V,p = PnVi/in + PpViflp (12) 

We also have 

dt5p + Vj{p5v') = Q (13) 

At this point we have two equations which are identical to the perturbation equations for a single fluid system. Of course, 
we are considering a two-fluid problem and there is a second dynamical degree of freedom. To describe this, it is natural to 
consider the difference in velocity. Thus, we introduce 

5w' = 5vl - Svi (14) 
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Combining the two Euler equations in the relevant way we have 

(1 - e)dt5w, + VrSp + 2B' eiju^^ Sw'" - 2BeijuCl^ e'^'^^iSwm = (15) 

We have defined 

5(3 = 5/ip — (5/in (16) 
which represents the (local) deviation from chemical equilibrium induced by the perturbations. We have also introduced 

£ = e„/a;p , B' ^ I - S'/xp , B = B/x^ (17) 

and we remind the reader that 

PpEp = Pn£n (18) 

Finally, we need a second "continuity" equation. To close the system, it seems natural to consider an equation for the proton 
fraction Xp. We then find that 

a...p + lv,[.p(i-.p)p..^]+.«^v,.p^o (19) 



2.3 Model equation of state 



Let us now consider the set of equations that we have written down. We see that the two degrees of freedom [5vi,5p\ and 
[5wi,5l3\ only couple "directly" through (19 1. For compressible models, the two degrees of freedom also couple "indirectly". 



since we can use the equation of state to relate [5p, 5(3] to [5p, 5xp\ 



Deciding to work with 5p and 5(3 [cf. [Andersson et al. (2008l], we use 

'dp 



and 



&Xp = 



dp 

dp 

dxp 
dp 



5p + 



8(3 



3(3 



5(3 



As a simple model, we shall consider an equation of state such that 

>^ =1 
dp) c1 

where Cs is the background sound speed. We combine this with a simple linear expression for the proton fraction; 



This leads to 



dxr 



a 

72 



Xp — ap 



and 



dXp 
~d(3 



2 2 

a p 



(20) 



(21) 



(22) 



(23) 



(24) 



= P 



dxp 
dp 



where we have used the relation ( Andersson fc Comer|[200T l 

.dp, 

to get 

■d£ 

In our model calculations we shall take ot — Q x 10~'^/pnuc where pn 

In a superfluid system the momentum of each component may not be parallel to it's velocity. Rather, it acquires a 
component along the relative velocity. This is evident from equation (|2|. This non-dissipative coupling is usually parametrised 



ap 



(25) 



(26) 



2.8 X lO'^* g/cm'' is the nuclear saturation density. 



in terms of the entrainment parameter Ex. One can show ( Prix et al.|2004 1 that this parameter is linked to the effective proton 
mass, rrip, according to 

For simplicity, we shall take Sp to be constant throughout the superfluid region of the star. Recent work suggests that, while 
this may not be a good approximation for an entire neutron star, it is approximately true if we consider a shell. From [Cliamel| 
(20081 we see that a reasonable range for the entrainment parameter is Sp « 0.2 — 0.8 



In order to compare our results to previous work, it is worth recalling that the entrainment parameter e used by |Lindblom| 
& Mendelll (|2000|) is related to £p by 
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Hence, their range of e ~ 0.02 — 0.06 corresponds, if we take a volume average of Xp, to £p « 0.45 — 0.85. 

To summarise; When we solve the r-mode problem, we will consider an equation of state represented by i) the overall 
density profile, represented by the sound speed c^, ii) the proton fraction Xp and iii) the entrainment parameter e-p. These 
quantities allow us to specify the background needed for the perturbation problem. This model is quite simplistic, and it 
would not be difficult to make it more realistic. However, our main interest is to explore how the r-mode results depend on 
the different parameters. This question is easier to address with a simple parametrised model. 



3 SLOW ROTATION ANALYSIS 



In order to determine the rotational corrections to the superfluid r-modes, we will apply the formalism of Saio 1 1982 1 to the 
problem. Thus, we consider corrections up to order 0(r2^) in the slow rotation approximation. We will also make the Cowling 
approximation, i.e. neglect perturbations of the gravitational potential, 5^. 

We start by considering a slowly, and uniformly, rotating star. It is well known that such a star will not be spherical and 
that a distorted potential surface can be written in the form 

r^a[l + e{a,e)] (29) 

Here e is a function of a and 6 which represents the deformation of the equilibrium structure from the background spherical 



state. The equilibrium physical quantities are functions only of a. Following Saio ( 1982 I; Smeyers & Martens ( 1983 1, we write 



the equations of motion in a frame corotating with the star, denoted by {g'}, starting from a static Cartesian frame, denoted 
by {x^}. Our coordinates in the rotating frame will be a spherical polar system, explicitly: 



X — a(l -I- e) sin ^ cos((j!) -I- Sit) 
x^ = a(l + e) sin 6 sm((j) + fit) 
x'^ = a(l + e) cos 9 

The metric in the new coordinates is then given by 



which, after linearizing with respect to e, leads to 



gab 



Sij 



dx' dx^ 



dab 



l + 2e + 2ai 




The connection coefficients, in the coordinate basis 



can be obtained from 



We shall, however, use the vector basis 



a^(l + 2e) 




d_ 

da ' 






a^(l-|-2e) sin^ 



+ 0(6^) 



ra 
be. 



1 ad ( I 

o5 \9bd,c + gad,b 



■ gbc,d) 



d_ 
da 



ld_ 
a 06 



for which the covector basis is 



= {da, ( 



i,asm 



(30) 
(31) 
(32) 

(33) 

(34) 

(35) 
(36) 
(37) 
(38) 



Note that the basis vectors are not orthogonal. 



Before we proceed, it is worth recalling that Smeyers & Martens ( 1983 1 have argued that the derivation of Saio's results 



is flawed, even though the final results for the rotational frequency correction are correct. Since this is an important point, 
we will outline how the second order slow-rotation perturbation equations should be derived following Saio's strategy. 

As we are considering linear perturbations on a rotationally distorted background, the first step is to calculate such 
a background, using the coordinates (o, 9, (j>) defined above. The equations that govern the background equilibrium are 
particularly simple (if we assume that neutrons and protons do not move relative to each other) 

1 dPja) ^ d<S>R{a) 
p{a) da da 

where <1>_r, defined in equation Q, includes the centrifugal terms. This is due to the fact that the coordinate a has been defined 
in such a way as to label the equipotential surfaces of "I>_h, which are also isobaric (and isopycnic as the background equation of 
state is barotropic) surfaces of the star, thus leading to all the background quantities being a function of a only [for a detailed 
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description see Tassoul ( 1978 1] . The Euler equations are thus identical to those one would obtain for a spherical background 
and the solution can be formally obtained by replacing the radial coordinate r with the variable a in the spherically symmetric 
solution. It is important to remember, however, that we are calculating physical quantities at a point P in the deformed star 
labeled by the coordinates {a,6,(ji) and not at a point Po in the spherical star labeled by coordinates (r = a,9,(j>). Note, in 
fact, that as the geometry is not spherical the measures of distance and volume change. For example, the mass of the star is 
now obtained from 



M : 



Je=o J 4,^0 



p{a)a sinl 



1 + 3e{a,e) + a 



de{a,9) 
da 



da dO d4) = Mo + SMn 



(40) 



where Mo is the mass of the spherically symmetric star and SMn is the difference in mass due to rotation. 



3.1 Perturbations 

Let us now consider linear perturbations on our rotationally deformed background. We shall express the perturbations in 
terms of a total displacement vector such that Svi = icj^f , and a difference displacement vector such that Swi = . 
We can then write equation (111 as: 



a^^t~-y'i5P^^5j-^+iaCt 







(41) 



p da 

where the 5 that denotes the Eulerian perturbations should not be confused with the Kronecker delta Sia- We are using iaC, 
to represent the Coriolis force, and the vector has the following components: 



Ca = 20 1 + 2e + a 



C2 



de 

'■IT 
oa 



smOQ 



Cl 



C±sin6» ( l + 2e + ( l + 2e + tan 



2f7 ( 1 + 2e + tan 6*^ ) cos 61^^ 



20, 



de 



The differential operator is the same as for spherical polar coordinates, i.e. 



da 



de 



a sine d(j) 



The continuity equation (fTsll becomes 



9e\ 



5p + (per) + pC+V^ [ie + a— J = 
We shall also need the linearized equation of state. In our model case, we have 



(42) 
(43) 
(44) 

(45) 
(46) 
(47) 



That is, in comparing our equations to Saio's results one should only consider the barotropic limit. 

Let us now consider the equations that govern the second degree of freedom, the difference in velocity. Equation ( 15 1 
takes the form 



where d has the following components: 



Ga = n 



1 + 2e + 2a 



da 



/i2 

sme 4_ 



+ ( 1 + 2e + + tanS^ ) cos6lsinef? 
da de ) 

Gg = n ( l + 2e + 2tan6l|^^ cosS^e* 

de de „de\ „ . 

1 + 2e + a— + cote— cosSsmef" 

da da do J 



G4, = 0(1 + 26) C 



Finally, we find that the second continuity equation ( 19 1 takes the form 

Sxp + -V'o [xp{l - Xp)p5~] + Xp(l - Xp)^~Vo ( 3e + 



da 



(48) 



(49) 

(50) 
(51) 

(52) 
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4 THE R-MODE PROBLEM 



We are mainly interested in understanding the efTect that mutual friction has on the r-mode instability. It is well-known that 
the r-modes have zero frequency and are purely axial in the non-rotating limit. Rotation breaks this degeneracy and leads 
to modes that are purely axial to leading order but which have a poloidal component of order 0{fl^). These modes have 
frequency 



(53) 



The r-mode frequency can be calculated, to second order in rotation, for a single fluid star using a number of different 
formalisms ( Andersson et al.|1999a [Lindblom et al.|1999 Yoshida fc Lee|2000 \. The superfluid problem has been solved both 



to linear (Andersson & Comer||2001| 



Prix et al. 2004 



Andersson et al.||2008[ ) and second ( [Lindblom fc Mendell||2000| |Lee fc 



|Yoshida|2003 1 order in rotation. To leading order one finds that the ordinary r-mode has frequency 

where, in fact, only the I = m modes exist. In addition, a constant density model supports a set of purely axial counter-moving 
modes with frequency ( Andersson et al.|2008 l 



2mn 



;(; + f)(f-e) 



B 



(55) 



Since the frequency in ( |55[ ) contains the term (f — e), it is easy to see that it can only be the solution for a global mode if 
e — Sp/ (1 — Xp) is constant. In the present analysis we will assume that Sp is constant. Hence, the second class of r-modes can 
only exist if we also assume that the proton fraction Xp is constant. If a::p is constant, the co-rotating and counter-rotating 
degrees of freedom are, in fact, completely decoupled also at second order in rotation. The r-mode is then exactly the same 
as for a barotropic single fluid star. In particular, r-mode solutions exist only for I = m. However, we are not generally 
considering a constant Xp (in fact we are assuming that the proton fraction scales linearly with the total mass density, which 
is not constant). Hence, we see there will be no counter-moving r-modes in our model. In fact, if Xp is not a constant, then 
the coupling between the degrees of freedom leads to the counter-moving r-mode becoming a general inertial mode, with a 



mixed toroidal/poloidal velocity field to leading order. These modes have been determined numerically by Lee & Yoshida 



(20031. If we consider a generic profile for Xp and work at second order in rotation the leading order r-mode will drive the 
counter-moving degrees of freedom, leading to mutual friction damping. This is the effect that we are interested in. 



4.1 Perturbation equations 



So far, we have written down the equations for a general perturbation problem on the rotationally deformed background. We 
will now focus on modes that are purely toroidal to leading order. That is, we write the total displacement vector as 



a 



Kim d 

sin d(j) ' 



d_ 
'89 



d z,, d 



Analogously, the difference displacement vector can be written 



a \ ' sin 9 d(f> ' 



, —ki. 



d 



+ E 



sin 9 d(j) , 



(56) 



(57) 



Note that we use uppercase variables for the co-moving degree of freedom and the corresponding lowercase variable for the 
counter-moving degree of freedom. 



As noted by Smeyers fc Martens (19831, the above form for the displacement vectors does not correspond to the usual 



decomposition in spheroidal and toroidal components. At second order in rotation, the above basis differs from the standard 
basis on a spherical background as the vector is no longer orthogonal to eg. This is, however, not a problem as we do 
not need to explicitly identify the poloidal and toroidal parts of the mode at second order in rotation. Such an identification 
would not be very useful anyway since the components are coupled at 0{Q^). 

We want to study the classical r-mode, i.e. a mode that to first order in rotation involves only the co-moving degree of 
freedom and which is purely toroidal in the slow-rotation limit. This identification is unique since our decomposition coincides 
with the standard one into poloidal and toroidal components at 0{i}). The requirement is that the leading term Kim is of 
order unity while the amplitudes of the "spheroidal" components are of order il^ for the total displacement. Meanwhile, 
following the first order results of Andersson et al. (20081, all components of the difference displacement must be of order Q^. 



As discussed in Appendix |Aj this ordering for the displacement vector, together with the frequency from equation (53 1, leads 
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to the following equations for the I = m r-modes; 



da 

da 
dsi+i 
da 
__^ dri+i 
da 



0+1 



^ - 3^ Si+i - ^0+1 + + + 2)Zi+i + 3imQi+iK,r^ ( 3£>2 + 
V 

(1 — U)C,l + l — Y^P'^' + ^ ^ 2iClUJLulQl+lKlm 



dlh 

da 



= {X — 3) si+i — Q+i 



V 



+ + + 2)zi+-, - CSi+i + 3imQi+iki^ ( 3D2 + 



' Xp{i -Xp) r 

(1 - U)ti+i - 2clujCjQi+iki^ [ilB - - 2ciujCjzi+i [mB + iB {{I + l)Qf+2 - {I + 2)(5f+i)] 



— —2iLuCjci 



l + l 



Qt+iKi, 



(58) 
(59) 
(60) 

(61) 
(62) 



-2iB Luujcilil + 2)Qi+iki,n + 2m(/ + 2)BujCjcxQi+\kir. 



+c^Lozi+^ |a;(/ + !)(/ + 2)(1 - e) - 2mCoB - 2iCoB [m^ + Q?+2(/ + + 4) + Q^i+^{1 + 2){l - 1)] | 
-2cia;iis,+i {mB - [Qf+2(/ + 4) ^ Qf+i(/ - 1) - l] } 



where 

klm 



■ ■ 

mrj 



IQi+i [Si+i + {l + 2)Zi+i] 
(mB - iZ^') Qi+i [s,+i + {l + 2)zi+i] 



In the last two expressions we have used 

2uj 



e = [lj ~ ujo) + 3D2 
Tj = {1 - e){ujo - ioi) 



l{l + l) 



IQi+i] - ^0 [5Q?+i - 1] 



(63) 
(64) 

(65) 
(66) 

(67) 
(68) 



It is easy to verify that these equations coincide with the results of Saio ( 1982 1 in the barotropic limit. In fact, in order to 
facilitate this comparison we have used essentially the same notation as Saio. Thus, we have defined the normalised mode 
frequency 



(69) 



and the following background quantities; 

47r / pa^ da 

Jo 
d log Mr 



Mr 
U 

V 

ci = ( 



dloga 

dlogP _ gap 
d log a P 
/a\3 M 
\r) M'. 

.3 \ 1/2 

0. ' 



X 

c 



\GM J 
dlog [pxp{l - Xp) 



dloga 
a dxp 
Xp{l — Xp) da 
e = Di(a) + D2(a)P2{cose) 

where P2{cos9) is the Legendre polynomial: 
We are also using the definition 



P2(cos6l) = (Scos^e- l)/2 



{I -m){l + m) 
{21 + 1){21 - 1) 



1/2 



(70) 
(71) 
(72) 
(73) 
(74) 
(75) 

(76) 
(77) 

(78) 
(79) 
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For later convenience, it is worth noting that this means that Qm = 0. 



4.2 Rotating n = I polytropes 

We will now restrict ourselves to the case where the density profile of the equilibrium configuration is that of an n = 1 
polytrope. This greatly simplifies the analysis, as the background quantities we are interested in can be obtained in closed 
form. Explicitly, we have 



p = Kp 



with 



K = 



2GR^ 



Introducing the dimensionless variable y = na/ R we then find that 

TvMo sin y 



y 



from which we can define 



Mr = [smy - 

TV 



y cos y] 



(80) 



(81) 



(82) 



(83) 



Furthermore, we can calculate, following the classical work of Chandrasekhar (19331, the rotationally induced deformation 

-2 



Di = DiQ' 



2 Moaipi ,2 

7r2 RMr ^ 



and 



P, n ~2 1 Moaip2 -.2 
D2 = D2u ="9^^- 



where 






tpi = 




sin y 

y 


1 - 




15 










y 





sm y cos y 

y 



(84) 

(85) 

(86) 
(87) 



The mass of the star is thus, from equation (401, given by 

Mo 



M ■ 



1 + 



In the following, when we quote the mass of the star we shall in fact be refering to the mass of the spherical star, Mq. Formally 
what we are doing is thus considering a sequence of stellar models with the same central density, but with a mass that varies 
with the rotation rate. The difference in mass does not, however, enter the r-mode frequency correction to order 0(f2^) and 
is thus irrelevant for the present discussion. 

Finally, we will also need the sound speed, which simply follows from 



= 2Kp 



(89) 



4.3 Decoupling the degrees of freedom 



Let us now examine the consequences of the ordering we assumed at the beginning of the analysis. From equation (62 I we see 
that if Kim is of order unity then (^im. (i.e. the pressure perturbation) must be of order 
that, as all components of the difference displacement are of order the variable r;„ 



Similarly, from equation ( 63 
(i.e. Sf3) must be of order fi* 



that in equation ( 59 1 the term involving t;. 



It follows 

is of higher order than the others and can be neglected. The consequence of this 
is that the equations for the co-moving degree of freedom are completely decoupled. Hence, they can be solved independently 
and then used as source terms for the counter-moving degrees of freedom. This approach has obvious advantages compared 
to previous fully numerical calculations ( Lindblom fc Mendell|2000 Lee fc Yoshida|2003 1. After all, we can now solve for the 
r-mode throughout the star, imposing regularity at the centre of the star and the vanishing of the Lagrangian perturbation 
of the pressure, Ap, at the surface. As a second independent step, we solve the equations for the counter-moving degree of 
freedom. Adopting this strategy it becomes straightforward to account for, for example, the temperature dependence of the 
superfluid gap and the associated variation of the superfluid region in the star. 

The solution for the comoving degree of freedom is particularly simple. Not only is it decoupled from the counter-moving 
degree of freedom, the equation for ("jm also decouples and takes the form: 



= (/ + 2-(/)0+i 

da 



(90) 
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leading to the solution 



0+1 _ ^ 



LOij I Mo / a\ '+1 



1 MrR \r) 



B 



1^ 

Mr 



One can use this solution to determine Si+\ from equation (58 1. This leads to 

+ + 4 - V,)Si+^ = -{Vg + /i)0+i 
where Vg — V/T and (for the I — m modes we are considering) 



h - 



1 M(a) fR 



uj^ M \ a 



{l + l){2l + 3)e 
KuiQ 



+ 3 3D2 + a 



dlh\ 
da J 



(91) 



(92) 



(93) 



with e defined in equation (A28I. To solve equation (921 let us first of all consider the solution to the homogeneous problem. 

C 



This is readily obtained and takes the form 



Si+i = 



pi/r„i+4 



(94) 



Note that this solution diverges both at the centre and the surface of the star. We next need to determine a particular solution 
to the problem. To do this let us write h as 

'R\ M{a) 



h ■- 



a J M 



-{he + ha) 



where 



ha 



(; + l)(2Z + 3)aj2 



(21 + 3)£>2 + a 



da 



Inspired by the solution to the homogeneous problem we make the following ansatz for the particular solution to ( 92 1 ; 

/(«) 



Si, 



pi/r„!+4 



This leads to 



The first term is thus of the form 



da 



= -(\4 + /i)//^a^+\,+i 



.,21+4, 



Vgp'/'^a-+\,+, = BG^ 
which is easily integrated. For the second part we need to be able to integrate, for the term proportional to /ic, 



hcp ' a Q+i = BV Khc — pa 



so that, for an n = 1 polytrope we require 



r 2i+2 , TtM ^i? 

/ px dx 
Jo 



21 + 3 ^na/R 



4i?3 



if) L 



21+1 . , nM fR 
y smydy= — 



21+3 



My) 



For the remaining part, involving ha, we require the integral 



I2 — S pa 



21+2 



{21 + 3)D2 + a 



dlh 
da 



which, using the explicit form for D2 in equation (851, gives 

5MR^ 



4tj21+2 



[fi{y) + f2{y)] 



where 

My) = 
My) = 



y^' siny 
sin y — y cos y 

ry 



[{3 - y'^ ) sin y-3y cosy] 



ry 

I X ~ [{3 — X ) sin X — 3a; cos a;] dx 
Jo 



(95) 

(96) 
(97) 

(98) 
(99) 
(100) 

(101) 
(102) 

(103) 

(104) 

(105) 
(106) 



This analysis is perhaps a little bit messy, but the result is very useful. Collecting the various results we can write the 
final solution to equation (|99l) as 



f{y) = 



BG fR 



21+5 



y 



TT^ (/ + l)^(2Z + 3) 



2(2; + 5) 16 



^2^1 {y) 



57r^(/ + l)^ 
8 



[My) + My)] 



(107) 
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I 




0"2 


1 


1.000 


0.3033 


2 


0.666 


0.4525 


3 


0.500 


0.4429 


4 


0.400 


0.4051 


5 


0.333 


0.3659 


6 


0.286 


0.3310 


7 


0.250 


0.3010 


8 


0.222 


0.2754 


9 


0.200 


0.2535 


10 


0.182 


0.2346 



Table 1. r-mode frequency corrections at second order in the slow-rotation approximation. The quantity (72 is defined by the relation 
a = <TQfl + cr2fl^ and, for the I = m case that we are considering, cro = 2n/(/ + 1). These results are in good agreement with the results 
of Lindblom et al. (1999), within a few percent for Z > 2, the greatest difference is approximately 10% for the I = m = 2 mode. These 
differences are most likely due to our use of the Cowling approximation. In fact for the I = m = 2 mode our frequency agrees to within a 
few percent with that obtained with the numerical code of Passamonti et al. (2008), in which the Cowling approximation is made, and 
to better than 1% with the results of Andersson et al. (2001). 



That is, we have an analytic solution to the problem. 

As we now have the full solution to the problem we can impose boundary conditions. First of all we require the solution 
to be regular at the centre of the star. This determines the constant C = in the homogeneous solution. The remaining 
condition is that the Lagrangian variation of the pressure vanish at the surface of the star, which in terms of our variables 
means that 



Ap = pga'^lCi+i - Si+i]Yi+i 



Given that our equation of state is such that p - 
that we must have 



as a^O (108) 
at the surface, we simply require (i+i and Si+i to be regular. This means 



fin) = 

Conveniently, this leads to an algebraic formula for the rotational frequency correction (T2, 



16 



I 



0-2 



{1 + 1)^(21 + 3) 



^21+5 



57r2(Z + 1)2 



/2(^) 



1 

l2 



(109) 



(110) 



2(21 + 5) 8 

The determined frequency corrections for the first few values of I are listed in Table [T] These results are in good agreement 
with the numerical results of [Lindblom et ar| ( |1999[ l, within a few percent for I > 2. The greatest difference is approximately 
10% for the I = m — 2 mode. The discrepancy should be entirely due to our use of the Cowling approximation. In fact for the 
I = m = 2 mode our frequency agrees to within a few percent with that obtained with the numerical code of |Passamonti et"^ 



(20081. Or results also agree well, with errors smaller than 1%, with the previous Cowling approximation result of 



Andersson 



& Kokkotas (20011 



4.4 The superfluid degree of freedom 

At this point we have reached two interesting conclusions. First of all, we have found that the rotational corrections for 
the r-mode can be obtained in closed form for n = 1 polytropes (in the Cowling approximation). This analytic solution 
will undoubtedly help us understand the nature of the r-modes and the effect of different dissipation channels better. The 
second result concers the fact that the superfluid, counter-moving, degree of freedom remains uncoupled also at this order 
of approximation. This is obvious since the frequency correction to order il"^ could be calculated without determining the 
counter-moving eigenfunctions. These are, in fact, only needed if we want to calculate the frequency correction to order 
As we shall see in the following this is exactly the order at which the mutual friction damping corrections appear explicitly 
in the frequency. 

To complete the r-mode solution at the current level of approximation we need to use the analytic solution to the co- 



moving problem to provide the source terms in equations ( A5 1 and ( A6 1 for the counter-moving degrees of freedom. When 



we consider this problem it becomes apparent that the decoupling of the two degrees of freedom is advantageous. In fact, it 
is now straightforward to study neutron stars with different models of superfluidity. This is an important improvement on 
previous work in this area, where it has generally been assumed that the entire neutron star core is superfluid. Our models 
will obviously still not be truly realistic, but we can (at least) compare results for different predicted superfluid energy gaps 
and varying neutron star core temperatures. 

In order to determine to what extent the neutrons and protons are superfluid/superconducting at a given core temperature 



(assumed to be uniform), we will consider the phenomenological gap models discussed by Andersson et al. (20051. It is 
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Figure 1. These figures show the critical temperatures for the onset of superfluidity of protons and neutrons as functions of radius in a 
steUar model with M = IAMq and _R = 10 km. We consider two models representing "strong" (left) and "weak" (right) superfluidity, 
respectively. For strong superfluidity we use models h and e of |Andersson et al.| | |200"5l l, while the weak superfluidity case corresponds to 
their models f and 1. The shaded region indicates the crust, the dynamics of which is not accounted for in our r-mode calculation. 



straightforward to use these models since they are given by analytic "fits" to the original results. The only caveat is that the 
fits are not valid at very low temperatures. One should also be careful to use the analytic gap functions only in the density 
region where they are relevant. 

We will compare r-mode results for two different models, intended to represent the range of possibilities. These models 
correspond to models [f,l] and [h,e] of Andersson et al. ( 2005 1 and represent what we will from now on refer to as "strong" and 
"weak" superfluidity, respectively. For our polytropic neutron star model this leads to the critical temperatures (as functions 
of the radius) shown in Figure [l] The results correspond to a model with mass M = IAMq and radius _R = 10 km (in 
the non-rotating limit). We take this as our canonical model. The data in Figure [l] confirms the expectation that below a 
critical temperature the superfluid region will always be a shell, between to radii Rmin{T) and -Rmax(r), where T is the core 
temperature. In order to simplify the analysis, we will only consider the core superfluids. This makes sense since we have not 
accounted for the elastic properties of the crust in the first place. Hence, we consider only the neutron superfluid and 
the ^So proton superconductor. To model the fluid dynamics of this system, we make the usual assumption that protons and 
electrons are locked electromagnetically. This leads to the system exhibiting two-fluid dynamics only in the region where the 
neutrons are superfluid. When the neutrons are normal, scattering with the electrons will lock the neutron fluid to the charged 
components. Given these assumptions the two-fluid region is such that Rmin{T) corresponds to the inner T = Tc point for the 
neutron gap. Meanwhile, 7?inax(T) corresponds to either the outer T = Tc point or the crust-core transition (Re), whichever 
is the smallest. In our analysis we assume that the crust-core transition takes place at a density oi p — 1.6 x 10^* g/cm^ 
( Haensel|200l| ). For our canonical star this leads to Rc ~9.32 km. 

We need to impose boundary conditions on the two-fluid equations at J?min(r) and -Rmax(r). Following Lindblom & 



Mendell ( 2000 1 and Lee & Yoshida ( j2003j we take the pressure perturbation to be continuous at each interface. This ensures 



that it is consistent to use the analytic solution for the co-moving part of the solution throughout the star. With these interface 
conditions we require 

e^(i?min)=0 and C-(^?max)=0 (111) 

Naturally, the r-mode mutual friction damping time will now depend on the size of the superfluid shell, i.e. the temperature 
of the star. 



We also want to be able to make a direct comparison with the results of Lindblom & Mendell ( 2000 1 and Lee & Yoshida 



(20031. To do this we consider a model such that the entire core is superfluid, with a critical density ps = 2.8 x 10^* g/cm^, 
corresponding to a transition radius Rs ~ 8.85 km in our model. We assume that the pressure perturbation is continuous at 
this interface, i.e. impose 



C-(i?s) = o 

The remaining condition in this case is regularity at the centre of the star. 



(112) 



4.5 The mutual friction damping timescale 

We will use the standard energy integral approach to assess the relative importance of different dissipation mechanisms and 
their impact on the r-mode instability. As discussed by, for example Andersson et al. (20081, this analysis is based on a 
functional E — Ek + Ep, where Ek and Ep represent the "kinetic" and "potential" energy, respectively. In the non-dissipative 
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Figure 2. This figure shows the r-mode damping timescale, tq, for a range of values of the entrainment parameter, £p. In order to 
facilitate a direct comparison with the results of [Lindblom fc MendeU] | |2000[ l and |Lee fc Yoshida| ( [2003^ , we consider a model where the 
entire core is superfluid. The transition density is taken to be ps = 2.8 X 10^* g/ cm'^ . The stellar model has mass M = IAMq and radius 
R = 12.533 km. We also assume that the mutual friction is due to electron scattering on vortex array, which gives B ~ while B is 
obtained from equation ( |117| l. The absence of "resonances", associating particular values of the entrainment with very short damping 
times, is notable in our results. 



case, the total energy is conserved. For the co-moving r-mode the problem simplifies considerably. The potential energy is 
higher order in rotation than the kinetic energy, so all we need is 

1 



p [\Sv\^ + (1 - e)xp{l ~ Xp)\&w\^] dV 



As discussed by |Andersson et al.| ( |2008| , the damping timescale for a mode r — —l/Q{uj) follows from 

1 _ 1 /dE 
T ^ ^2E \~dt 



Focussing on the mutual friction dissipation, we find from equations (111 and ( 15 I that 

dE __ 

dt ^ 



(113) 



(114) 



(115) 



At the current level of approximation, we cannot determine the mutual friction damping of the classical r-mode directly 
from the imaginary part of the mode frequency. This would require calculating the frequency up to order 0{^l^), not 0{^l'^) 
as we have done here. This is easily seen if we note that, as the mode is purely co-moving to first order, Sv ~ 0(f2) and 
5w Hi O(n^). Thus, the energy scales as 



E ^ {SvY ^ 0{n^) while 



dE 
dt 



and, from equation (114 1, we see that Q(lo) ~ 0{D, ). This means that we, in fact, have to resort to the energy integral 



estimate. Since the evaluation, to leading order, of the integrals in (1131 and (1151 only requires the eigenfunctions up to 
order 0{n^) it can be performed within the present scheme. Due to the expected scaling we will focus on the quantity tq 
defined by the relation 



1 

To 



n 



(116) 



where p is the mean density of the star. 

Having established the strategy, let us first of all consider the model where the superfluid extends all the way from the 
centre of the star up to Rs- In order to compare our results directly to |Lindblom fc Mendell| ( |2000[ ) and |Lee fc Yoshida| ( |2003| ) 
we take the mass and radius of the star to be M = 1.4M0 and R = 12.533 km. We also assume that the main cause of mutual 
friction is electron scattering off the vortex array. This assumption puts us firmly in the "weak" drag regime where we can 
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Figure 3. This figure sliows the r-mode damping timescale, tq, calculated for ep=0.6 and 7?. = 1 for the strong superfluidity model. 
The results illustrate the importance of using the dissipative mode-solution in the evaluation of the energy integrals. The dashed curve 
shows the timescale obtained by integrating the undamped eigenfunctions. The other two curves show the effect of introducing first B' 
(dash-dot) and then also B (solid) in the calculation. 



neglect the B coefficient. Meanwhile, for B we have the result of 



Andersson et al. 



(20061; 



5 = 4 X 10" 



1/2 



7/6 



.0.05/ 



1014 g/cm-' 



1/6 



(117) 



Results for the damping timescale to due to mutual friction, with varying Sp, are presented in Figure [2] These results confirm 
the main conclusion of Lindblom & Mendell (20001 and Lee & Yoshida (20031. Mutual friction due to electron scattering off 
the vortex array is too weak to affect the r-mode instability significantly. There is, however, one important difference between 
our results and the previous studies. We do not find any value of £p for which the damping timescale becomes short, cf. the 
resonances discussed by [Lindblom fc Mendell ( 2000| and 



Lee fc Yoshida| ( |2003[ ). The absence of these resonances is most 



likely due to the fact that, following Andersson et al. (20081, we have imposed that the mode is purely axial to leading order. 



That is, we neglect higher order terms that would couple the counter-moving motion back to the co-moving motion. In effect, 
we have a priori ruled out the resonant behaviour found by Lindblom & Mendell (20001. Moreover, we cannot have avoided 
crossings with the superfluid inertial modes as discussed by Lee & Yoshida (20031. Once we assume that the r-mode is purely 
axial to leading order, our analysis becomes completely oblivious to the fact that other mode solutions to the general problem 
may exist. 



Another important difference between our analysis and the work of Lindblom & Mendell (20001 and Lee & Yoshida 



( 2003 1 is that we explicitly keep the mutual friction terms in the equations of motion. The damping timescale we calculate 
is obtained by using the full 0(f2^) solution for the problem in the integral in equation (114 1, rather than integrating the 
solutions to the inviscid problem. In the weak drag limit, e.g. when we consider (1171, this does not affect the results at all. 
However, by keeping the mutual friction force in the perturbed equations of motion we are no longer restricted to the weak 
drag limit. We can consider the entire range < 7?. < oo. The capacity to consider the strong drag regime may, in fact, be 
quite important. Dynamics in the strong drag regime has only recently been considered, see e.g. Glampedakis et al. (2008a I, 
and the results show that the problem has interesting features. Hence, we will consider the r-mode problem outside the range 
TZ <^ 1. Motivation for studying the problem for TZ^ 1 comes from the possibility that the interaction between fluxtubes and 
neutron vortices may be efficient ( [Ruderman et aL]|1998[ |Link||2003 20061. The problem would also be in this strong drag 
regime if a fiuxtube cluster is associated with each neutron vortex ( Sedrakian fc Sedrakian|1995 I or in the presence of strong 
vortex pinning ( Shaham|1977 1. 

Because of the "symmetric" dependence of the mutual friction coefficient B on the drag TZ, one may expect the damping 
timescales to be similar for values TZ and 1/7?.. The damping should be most effective for TZ ^ 1. The results in Figure [3] 
illustrate the importance of keeping the mutual friction force terms in the perturbed equations of motion in this efficient 
friction regime. In the figure we compare the mutual friction damping timescale determined from the energy integral using 
eigenfunctions where we either keep both the B and the B' terms or (artifically) set one or both of them to zero. The evidence 
is clear. For TZ ~ 1 the inviscid solution does no longer lead to a good approximation. 
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5 THE R-MODE INSTABILITY WINDOW FOR STRONG MUTUAL FRICTION 



It is interesting to consider the case TZ ~ 1 further, even though it is somewhat extreme (and possibly only of academic interest). 
After all, it leads to the strongest possible mutual friction damping. Moreover, because it corresponds to Z3 ~ Z3' ~ 1/2 it 
provides insight into the dynamics of problems where the two terms in the mutual friction force are of similar importance. 

Let us consider the effect that a strong mutual friction would have on the r-mode instability. To do this we need to 
account for the various viscous processes which, at a given temperature, act to damp the gravitational-wave driven instability. 
First of all it is easy to see that the r-modes will only be unstable in a certain temperature range. At temperatures below 
T ~ 10^ K shear viscosity will always suppress the instability, while for temperatures above T ~ 10^^° K (at which no part of 
the star is expected to be superfluid) bulk viscosity will prevent the mode from becoming unstable. Our aim is to establish to 
what extent the mutual friction can further restrict the range in which the mode is unstable. 

The r-mode instability window is usually illustrated by the critical rotation period above which the mode is unstable as 
a function of temperature. The relevant critical rotation rate is obtained by solving for the roots of 

1111 

\ h ■ 

where Tgw is the growth timescale for the instability due to gravitational-wave emission. For an n 
approximated by ( Andersson &: Kokkotas|2001 1 

M \"V R \~^' f P 



+ — =0 

Tmf 



(118) 

1 polytrope it is 



Even though this estimate was obtained for a single fluid star, it will remain a good approximation for the ordinary r- 
modes of a superfluid star. This is obvious since the gravitational-wave emission is entirely due to the co-moving degree of 
freedom ( Andersson et al.|2008 1. Meanwhile, Tbv represents the damping time due to bulk viscosity. The bulk viscosity will be 
suppressed in superfluid regions. Hence, it is only active in regions where the star is not superfluid. In general, this leads to 
a complex problem. However, in the case of npe-matter where only the modified URCA process is relevant the bulk viscosity 
damping is not important at low temperatures. Hence, we essentially only need the bulk viscosity at temperatures above the 
superfluid transition. That is, we use ( Andersson fc Kokkotas|2001 1 



Tbv 



2.7 X 10' 



M 



R 



^1.4M0 7 VlO kmy \lmsj \lO^K J " '•^^^^ 
It should be noted that the role of bulk viscosity will be different for hyperons and deconfined quarks because of a different 
scaling with temperature, see e.g. Nayyar & Owen ( 2006[ ) for a recent discussion. 

The damping timescale, Tmf , for mutual friction damping, and r^v, the timescale for damping due to shear viscosity require 
a more detailed discussion. In the case of mutual friction we are interested in the temperature dependence of the damping 
timescale. To investigate this we consider the realistic gap models described in the previous section which, as they predict 
how the extent of the superfluid region varies with temperature, allow us to predict the temperature dependence of Tmf itself. 
Meanwhile, the shear viscosity damping also depends on which layers of the star are superfluid, as the dominant effect will 
be due to different processes in superfluid and normal fluid regions. In the region where the neutrons are not superfluid (the 
normal fluid region) we assume that the main contribution is due to neutron-neutron scattering. This leads to a viscosity 
coefficient ( Andersson fc Kokkotas|2001 1 

9/4 / m \ -2 



r;nn = 2 X 10" 



g/cm s 



(121) 



IQiSg/cmV VlO^'K^ 

If only the neutrons are superfluid then the dominant process will be electron-proton scattering which leads to the coefflcient 



(Andersson et al. (2005l) 



rjcp = 1.8 X 10' 



13/6 / 



13/6 



/TV 



g/cm s 



(122) 



0.017 V10"g/cmV \10^KJ 
If, on the other hand, both neutrons and protons are superfluid the dominant process is electron-electron scattering ( Andersson] 
et al.|2005 l, which leads to the coefficient ( Andersson fc Kokkotas|2001 1 



7?co = 6 X 10' 



g/cm s 



(123) 



VlOiSg/cmsy VlO^K, 

Once we have determined the dominant viscosity agent in each region of the star we can compute the damping timescale from 
the energy integral 

'dE' 



dt 



= -2 



.dV 



where the shear Saab is defined as 



\/b^t - 2gat^,e+) 



(124) 
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Figure 4. The r-mode instability window for ep = 0.6, calculated as a function of core temperature T and for 7^ = 1 in the "weak" 
(left) and "strong" (right) superfluidity cases. We consider a star with M = IAMq and R = 10 km. The dotted line indicates the shape 
the instability window has if we ignore mutual friction, while the dashed line indicates the effect that the Ekman layer at the base of 
the crust would have on the instability region. The results show that, when 7?. = 1, the mutual friction is strong enough to suppress the 
r-mode instability completely as soon as the core becomes superfluid. 



The obtained results for the instability window are shown in figures [4] [5] and [6] for a range of drag parameters in the 
7?. ~ 1 regime, and the "strong" and "weak" superfluidity models. The data corresponds to a A'l = IAMq and R — 10 km 
neutron star. In the figures, we express the critical angular velocity in terms of ilo = ypJ^Gp where p is the average density, 
and assume that the star cannot spin faster than the breakup limit, taken to be f2b ~ 2/3Slo. 

Note that, as we are considering a purely fluid star and do not account explicitly for the elastic crust we have not included 
the damping timescale due to the Ekman layer which is expected to form at the crust-core boundary. To indicate the effect 
that the Ekman layer may have on the instability window, we use the rough estimate 

p \ 1/2 



3 X 10' 



VlOSKy/ Vims 



(126) 



We arrive at this estimate by taking the simple constant density estimate of Andersson & Kokkotas (2001 1 for a M = IAMq 
and i? = 10 km neutron star, corrected for a "slippage" factor 5c=0.05, as defined by Glampedakis & Andersson (2006b I. 
In fact, it has been shown by Glampedakis & Andersson (2006a I that one should expect the constant density estimate to 
only differ by factors of a few from the result for a stratified model. Hence, it should be a reasonable approximation for our 
discussion. 

Perhaps not very surprisingly, our analysis demonstrates that the mutual friction damping can have a significant effect 
on the r-mode instability window in the extreme case 7?. ~ 1. It is, however, interesting to note the impact on the instability 
window of the two superfiuidity models. The results in figures [5] and [6] show that, for values of the drag parameter in the 
range 0.005 < TZ < 500, mutual friction may have a significant effect on the r-mode instability. Comparing the two figures, 
we see that the "strong" or "weak" superfluidity models lead to significantly different results in the strong friction regime 
(keeping the mass and radius of the star fixed). As the star cools below the superfluid transition temperature the mutual 
friction timescale initially decreases (the critical frequency increases). However, at lower temperatures the damping time may 
increases again (around T ~ IQ^K). This would lead to a second minimum in the instability curve. For parameter values 
outside the range 0.005 <TZ < 500, the mutual friction would not be the dominant dissipation mechanism for the superfluid 
r-modes. 

It should be stressed that we have focussed on cases where both B and B are relevant. In the very strong drag regime 
(i.e. TZ^ 1) one will have B ~ 1 but B <C 1. In this limit, one again finds that the mutual friction damping of the r-modes 
is irrelevant. 



6 CONCLUDING REMARKS 

In this paper we have re-examined the problem of mutual friction damping in rotating superfiuid neutron stars and its effect 
on the r-mode instability. Our analysis differed from previous efforts in that we did not a priori assume that the drag on the 
superfluid vortices is weak. By including the mutual friction force in the equations of motion, we were able to (for the first 
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Figure 5. The r-mode instability window for Sp = 0.6, calculated as a function of core temperature T and for a range of drag parameters 
TZ. We consider a star with M = 1.4Mq and ij = 10 km and "strong" superfluidity. The dotted line indicates the shape the instability 
window has if we ignore mutual friction, while the dashed line indicates the effect that the Ekman layer at the base of the crust would 
have on the instability region. The results show that, for values of the drag parameter in the range 0.005 < 7?. < 500, mutual friction 
may have a significant effect on the r-mode instability. It is particularly interesting to note the local minimum that may be present at 
low temperatures. This feature arises due to the effect shown in figure [s] 



time) consider r-modes in the strong drag regime. We calculated these modes to second order in rotation, thus extending the 



first order results of Andersson et al. ( 2008 1 



Our mode analysis focussed on solutions that are purely axial to leading order. For the classical r-mode, which is expected 
to lead to the fastest growing gravitational-wave instability, this assumption has the advantage that the equations for the 
co-moving degree of freedom decouple and serve as source for the counter-moving motion in the superfluid region of the star. 
For the particular case of an n = 1 polytrope, and in the Cowling approximation, we determined a useful analytic solution for 
the co-moving motion. The existence of this solution, which is relevant also for the single fluid barotropic r-mode, simplifled 
the solution of the superfluid problem considerably. 

In the superfluid problem there may also exist a class of counter-moving r-modes. These modes are, however, in general 
not purely axial to leading order. For a stratified neutron star model, the velocity field of these modes acquire a leading 
order polar component. Such modes cannot be determined within our current framework. This means that we cannot confirm 
the proposed existence of "avoided crossings" between the two classes of modes. Such crossings have been suggested as the 
explanation for the sharp resonances with very short mutual friction damping time scales found by |Lindblom fc Mendell| 
(20001 and Lee & Yoshida (20031. Our results do not show such resonances. Whether this is an artefact of our approximation 
scheme is not clear. We see no evidence that our approximation fails for particular values of the entrainment, as one might 
expect if the r-modes change character as the entrainment is varied. This issue requires further attention. 

We have, however, calculated the counter-moving r-mode in the only case where it can exist, when the star is not stratified, 
see Appendix B. In this case it is possible to obtain, once again, an analytic solution. This solution, which was presented here 
for the first time, has already been used as a test bench for numerical simulations of superfiuid neutron stars ( [Passamonti et| 
al.|2008| . 



For the classical (mainly co-moving) r-mode, our results confirm that mutual friction due to the most commonly considered 
mechanism (electron scattering off the vortex array) acts on a timescale that is too long to suppress the r-mode instability. 
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We conclude that, in the weak drag regime, the mutual friction is not the leading damping mechanism for these modes. This 
agrees with the conclusions of Lindblom & Mendell (20001 and Lee & Yoshida (20031. 

We have also considered, for the first time, the effect of the mutual friction on r-modes in the strong drag regime. In 
the extreme limit, this leads to the expected result that the damping timescale is (again) very long. Hence, mutual friction 
does not damp r-modes effectively in the TZ^ 1 regime. In order to understand the effect that a strong mutual friction may 
have, we have also studied the intermediate regime where TZ ~ 1, when S ~ S ~ f/2. This regime could be relevant if there 



exists regions in the star where the interaction between neutron vortices and fluxtubes is efficient (Ruderman et al. 1998 



pnk|[2003l[2006| . In this case we find that mutual friction can have a distinct effect on the r-mode instability window. The 
instability may, in fact, be completely suppressed below the superfluid transition temperature. Moreover, the mutual friction 
damping depends strongly on the model for the superfluid energy gap, which determines the critical temperature below which 
neutrons or protons become superfluid. We have compared results for two typical models, taken from Andersson et al. (20051, 
representing "strong" and "weak" superfluidity, respectively. This analysis lays the foundation for studies of realistic neutron 
star models, where the size of the superfluid regions depends explicitly on the stars temperature. 

With this analysis we have prepared the ground for a more detailed study of neutron stars with exotic cores, e.g. 
dominated by hyperons or deconfined quarks. In each of these cases, one would expect superfluidity to be relevant. In the 
case of hyperons, the problem is likely to require additional "fluid" degrees of freedom. Naively, one would expect the neutron 
and proton superfluid/superconducting mixture to coexist with a neutral A superfluid and a E~ superconductor. Assuming 
electromagnetic coupling between the charged components, this would leave three distinct hydrodynamical degrees of freedom 
(at zero temperature). Neutron stars with deconfined quark cores may also require additional degrees of freedom. Perhaps 
the simplest possibility corresponds to the so-called CFL phase ( Alford et al.||1999 l for which, at low temperatures, one has 



to consider a superfluid coupled to phonon excitations (Mannarelli et al 
that is formally equal to superfluid He** 



Andersson & Comer 



20081. This should lead to a dynamical problem 



(20081 for a recent discussion. More complicated phases, 
requiring the inclusion of Kaons, either as thermal excitations or as a condensate, will also need to be considered. It would be 
very interesting, and perhaps important, to understand better to what extent multifluid dynamics plays a role for such exotic 
systems. 
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APPENDIX A: DERIVING THE R-MODE EQUATIONS 



In this Appendix we provide the derivation of the equations that determine the rotational corrections to the r-modes. We take 
as our starting point the general slow-rotation perturbation equations from Section [3. 1| Focussing on modes that are purely 
toroidal to leading order, we assume that the total displacement vector takes the form 



= 



K,^ d 



d_ 
09 



d z,^ d 



'd9'- sin6, 



Meanwhile, the difference displacement vector can be written 



= 



sin 9 d(f> ' 



d 



d z. 



sin 9 



(Al) 



(A2) 



We use uppercase variables for the co-moving degree of freedom and the corresponding lowercase variable for the counter- 
moving degree of freedom. 

Our main interest is in the "classical" I = m r-mode, a mode that to first order in rotation involves only the co-moving 
degree of freedom and which is purely toroidal. This leads to the requirement that the leading term Kim is of order unity 
while the amplitude of the spheroidal components is of order Q,^ for the total displacement. Moreover, the first order results 

1— 1 I f n 

of Andersson et al. (2008 ) show that all components of the difference displacement must be of order . With this ordering 
for the displacement vector, and for the frequency given in equation (53 1, we have the equations 

dSi/n 



E 



da 



da 



E 



V 



V 



r 

+3imKim { 3D2 + a 



dih 
da 



cos^F, 



= E 



(1 U^C,l/^ ^ ^p'T~L'fJ 



— 2iciuJLjKim sin 



(A3) 
(A4) 



These relations follow from the equations for the total displacement, (111 and ( 13 1. From the corresponding equations for the 
difference, i.e., (|15[) and (|19[l, we get 



ds 



^ da ^ 



1 V 

{X - 3) S^^ - Ci^^j- r— + V{V + l)z„ij, - CS„fj 

(1 - xp) r 



(A5) 



da 



E 



(1 — (7)r„^l^j'— 2kimCiujuj ( iB sin 



-mB cos 9Yi^ 



2ciijjCjZvu. ( mB -f sin S cos 6* 



89 



+c^Los,^, {lo{1 - e)Y^ + 2iCoB{cos' 9 - l)Yj; 
In writing down these equations, we have defined the variables 

Sp 
P 



X^t Y" - - 



ga 



5(3 



(A6) 

(A7) 
(A8) 



The various background quantities are defined by equations (70|-(77l 

components of the co-moving I 

^ u{u + l)C,r.^,Y^ + 2iciUjCbKi 



Combining the 9 and (j) components of the co-moving Euler equation we also obtain two algebraic relations 

dYr 



l{l + l)cos9Yr +sine- 







(A9) 



and 



Kimuj \ l{l + 1)(1 + 2e){uj - iuo)Yr + 6D2 



to cos H sm f 



89 



2iLuCj J2 \ - '^{'^ + 1)^-/^1 cos 9YI' + {S^^ - Z,^) sin 



mCj{5 cos 9^ - 1)Y" 



.dY^ 



89 



(AID) 
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where 



2mCj 

uT+rj 



(All) 



Analogously, we obtain from the difference Euler equations 



1(1 + 1) cos evr + sine - 



de 



uju{u + 1) (1 - e) - 2£l}mS F^" - 2i(iS(mV„'' + + 1) cos^ OY^ + 2 sin 6 cos 



+ 2mBciLuujkim 



„2 /JV'A' 



2 COS 9Yi^ + sin 



-2cia;a) ^ s^^ - iB 



fly 

COS e sin + (3cos^e» - l)y^, 



(A12) 



and 



ki^ulil + 1)(1 + 2e)(l -e)(a;-LJi)y;™ = 2iuCj^ [(26 - imS)s^^ - + 1)2^,,] cos^yj^ 

r -' -' - 1 fly 

+ P ^ ~ imB)z^^j sin 6*-^ 



where 



/(Z + l)(l-e) 
We can now use the standard recurrence relations 
n9Yi 



^j— — r |2mS' + 2iB [1(1 + 1) - m^] j 



96' 



= «Qi+iiiTi-a + i)Qi^!"i 



cossyr = Q,+^Yi'];:, + QiY,^, 

with Qi defined in ( |79| |. Repeated use of these relations leads to 
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(A14) 

(A15) 
(A16) 
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[IQf+i -(1 + l)Qf]Yr + lQi+iQi+2Yt^2 
-(l + l)QiQi-iYi"^2 



(A17) 
(A18) 



cos" 9Yr = (Qf+i + Qt)Yr + Qi+iQi+2Y/^2 + QiQi-iY{I!2 
This way we arrive at the equations that that determine the next order slow-rotation correction to the I — m r-mode; 
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where 



K,^ = -^/(3,+i[5',+i + (; + 2)Z,+i] 
me 

0^0 



mrj 
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In the last two expressions we have used 
2lu 



= {U! — UJo) + 

= (1 — e)(a;o — t^i) 



/ + 1 



f+i - 1^0 [5Qf+i - 1] 



The above equations are identical (58l-(68l in Section 
slow-rotation. 



4.1 



(A28) 
(A29) 

They completely specify the r-mode problem at second order 



APPENDIX B: THE COUNTER-MOVING R-MODE 



In this appendix, we will examine the "superfluid" r-mode, a mode that to first order in rotation involves only the counter- 
moving degrees of freedom and which is purely toroidal. This leads to the requirement that the leading term kim is of order 
unity while the amplitudes of the spheroidal components are of order Q,^ for the difference displacement. Moreover, all the 
components of the total displacement are of order fl^. We have already mentioned that such a mode cannot exist unless the 
proton fraction is constant ( Andersson et aL||2008 1 . Nevertheless, these modes may be of some interest. Hence, it is worth 
determining the relevant slow-rotation corrections. 

Given the ordering for the displacement vector, together with the frequency from (53 1, we have the following equations 
(for the I — m mode) 
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We can now show that one cannot in general have a pure countermoving r-mode, as in this case one has kim of order 
unity and thus, from equation (B6I one finds that Tim is of order Q,^. This means that in equation (B2| the terms involving 
Om must be of the same order, but this would lead, from equation ( B5 1 to one or more components of the total displacement 
vector being of order unity. The mode would thus no longer be a pure counter-moving r-mode, but a more general inertial 
mode. However, this argument fails if we assume constant density profile. In this case there is no coupling term in equation 
(B2l and one can solve the problem by defining an enthalpy-type variable, see Andersson et al. (20081. 

Hence, we consider an equation of state for which the two fluids are decoupled. In particular, we shall consider the 
analytical equation of state of Prix et al. (20041; Passamonti et al. (20081. Then the energy functional takes the form 



where A^y is given by 
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and Xp, a and K are constants. In the background one has p ; 
one finds that 

&P = 2KpSp, 



Kp and a constant proton fraction, while for the perturbations 



2K{l + a) 



1- 



rpSXp . 



(B13) 



L,p - (1 + cr)a;p] 

We now see that the couphng term in equation ( |19[ | vanishes, as Xp is constant, and there is also no coupling from the equation 
of state as i) the pressure is a function of the density only, and ii) /9 is a function of Xp only. 

Let us consider the inviscid problem. This is of immediate interest since we can test the solution against the recent 
time-evolutions carried out by Passamonti et al. (20081. Neglecting the mutual friction, the equations for the mode we are 
interested in take the form 
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The equations for the co-moving degree of freedom are completely decoupled and take the same form as in equations ( |90[ l and 
(92 1. The frequency of the counter-moving r-mode can be written as 



(7* = (JoH + a|57 



with ( [Andersson et al.|2008 l 
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(BIT) 
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(1 - e) l{l + 1) 

In order to determine the correction (t| we adopt the same strategy as for the classical r-mode and use the solution to equation 
(B14I to find a particular solution to equation (B15|. Skipping the details of the calculation, which is essentially the same as 



in section |473| we obtain the relation 
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Here we have assumend that the entire star is superfluid, and imposed the boundary condition A/3 = at the surface, in 
order to be able to compare directly to the numerical results of Passamonti et al. (20081. This comparison leads to a very 
good agreement (see figure 2 of Passamonti et al. (2008l). 
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